http://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.fft.html#numpy.fft.fft
Compute the one-dimensional discrete Fourier Transform.
This function computes the one-dimensional n-point discrete Fourier Transform (DFT) with the efficient Fast Fourier Transform (FFT) algorithm [CT].
Parameters:
a : array_like
Input array, can be complex.
n : int, optional
Length of the transformed axis of the output. If n is smaller than the length of the input, the input is cropped. If it is larger, the input is padded with zeros. If n is not given, the length of the input along the axis specified by axis is used.
axis : int, optional
Axis over which to compute the FFT. If not given, the last axis is used.
Returns:
out : complex ndarray
The truncated or zero-padded input, transformed along the axis indicated by axis, or the last one if axis is not specified.
Raises:
IndexError :
if axes is larger than the last axis of a.
See also
numpy.fft
for definition of the DFT and conventions used.
ifft
The inverse of fft.
fft2
The two-dimensional FFT.
fftn
The n-dimensional FFT.
rfftn
The n-dimensional FFT of real input.
fftfreq
Frequency bins for given FFT parameters.
Notes
FFT (Fast Fourier Transform) refers to a way the discrete Fourier Transform (DFT) can be calculated efficiently, by using symmetries in the calculated terms. The symmetry is highest when n is a power of 2, and the transform is therefore most efficient for these sizes.
The DFT is defined, with the conventions used in this implementation, in the documentation for the numpy.fft mo
In [2]:
import numpy as np
import math
import iris
In [5]:
%matplotlib inline
In [12]:
experiment_ids = ['djzny']
pp_file_path='/nfs/a90/eepdw/Data/EMBRACE/'
diag = '408_on_p_levs'
p_lev=925.
In [13]:
for experiment_id in experiment_ids:
expmin1 = experiment_id[:-1]
fu = '%s%s/%s/%s_%s.pp' % (pp_file_path, expmin1, experiment_id, experiment_id, diag)
cube = iris.load_cube(fu)
cube_slice = cube.extract(iris.Constraint(pressure=p_lev))
In [14]:
cube.coord('pressure')
Out[14]:
In [15]:
cube_slice
Out[15]:
In [98]:
n=cube_slice.shape[0]*2
#int(math.ceil(np.sqrt(n))
In [134]:
def NextPower2(value):
powof2 = 1
# Double powof2 until >= val
while( powof2 < value): powof2 *= 2;
return powof2
In [135]:
powof2 = NextPower2(cube_slice.shape[0])
In [218]:
DiurnalFreq = 1/freq[np.where((1/freq<=28.) & (1/freq >=0.))]
DiurnalTransformAmps = np.sqrt(sp.real**2 + sp.imag**2)[np.where((1/freq<=28.) & (1/freq >=0.))]
AmpSort = np.argsort(DiurnalTransformAmps,axis=0)
DiurnalFreqSort = DiurnalFreq[AmpSort]
DiurnalAmpsSort = DiurnalTransformAmps[AmpSort]
In [220]:
DiurnalAmpsSort.shape
Out[220]:
In [215]:
DiurnalFreq[23]
Out[215]:
In [208]:
DiurnalTransformAmps[AmpSort]
Out[208]:
In [227]:
AmpSort[-1].shape
Out[227]:
In [209]:
DiurnalAmpsSort = DiurnalTransformAmps[AmpSort,:,:]
In [211]:
DiurnalAmpsSort.shape
Out[211]:
In [226]:
plt.contourf(cube_slice.coord('grid_longitude').points, cube_slice.coord('grid_latitude').points, DiurnalTransformAmps[3])
plt.colorbar()
plt.title(DiurnalFreq[3])
plt.show()
In [225]:
plt.contourf(cube_slice.coord('grid_longitude').points, cube_slice.coord('grid_latitude').points, DiurnalTransformAmps[2])
plt.colorbar()
plt.title(DiurnalFreq[2])
plt.show()
In [224]:
plt.contourf(cube_slice.coord('grid_longitude').points, cube_slice.coord('grid_latitude').points, DiurnalTransformAmps[23])
plt.colorbar()
plt.title(DiurnalFreq[23])
plt.show()
In [228]:
plt.contourf(cube_slice.coord('grid_longitude').points, cube_slice.coord('grid_latitude').points, DiurnalTransformAmps[AmpSort[-1]])
plt.colorbar()
plt.title(DiurnalFreq[23])
plt.show()
In [230]:
DiurnalTransformAmps[AmpSort[-1]].shape
Out[230]:
In [184]:
DiurnalTransformAmps[np.where((1/freq<=28.) & (1/freq >=0.))].shape
Out[184]:
In [147]:
freq = np.fft.fftfreq(powof2)
sp = np.fft.fft(cube_slice.data, n=powof2, axis=0)
half_transform = powof2/2
#half_transfrom = 246
#sort_by_amplitude = np.argsort(sp[:half_transform,0,0].real)
#1/freq[sort_by_amplitude]
In [172]:
sp.shape
Out[172]:
In [ ]:
sort_by_amplitude = np.argsort(sp[:half_transform,0,0].real)
In [171]:
plt.plot(1./freq[:half_transform], np.sqrt(sp[:half_transform,0,0].real**2 + sp[:half_transform,0,0].imag**2))
#plt.ylim(0,5000)
plt.xlim(0.,24.)
#plt.ylabel()
plt.xlabel('Frequency (hours)')
plt.show()
plt.plot(1./freq[:half_transform], np.sqrt(sp[:half_transform,0,0].real**2 + sp[:half_transform,0,0].imag**2))
#plt.ylim(0,1000)
#plt.xlim(0.,24)
#plt.ylabel()
plt.xlabel('Frequency (hours)')
plt.show()
plt.plot(sp[:,0,0].real, sp[:,0,0].imag)
#plt.ylim(0,1000)
plt.xlim(-1500.,1500.)
#plt.ylabel()
plt.xlabel('Frequency (hours)')
plt.show()
In [116]:
sort_by_amplitude
Out[116]:
In [121]:
freq
Out[121]:
In [122]:
len(sp[:,0,0].real)/2
Out[122]:
In [141]:
Out[141]: